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Abstract. The high density medium that characterizes the central regions of starburst galaxies and its power to accelerate 
particles up to relativistic energies make these objects good candidates as y-rays sources. In this paper, a self-consistent model 
of the multifrequency emission of the starburst galaxy NGC 253, from radio to y-rays, is presented. The model is in agreement 
with all current measurements and provides predictions for the high energy behavior of the NGC 253 central region. Prospects 
for observations with the HESS array (and comparison with their recently obtained data) and GLAST satellite are especially 
discussed. 



1. Introduction 

Starburst galaxies (and in general star forming regions) are ex- 
pected to be detected as y-ray sources. They have both a large 
amount of target material and, due to the presence of super- 
nova remnants and the powerful stellar winds of their numerous 
young stars, multiple shocks where to accelerate particles up to 
relativistic energies. Pion decay production of y-rays, usually 
dominant under such conditions, is thought to produce signifi- 
cant y-ray fluxes. 

Approximately 90% of the high-energy y-ray luminosity 
of the Milky Way (-1.3 xlO^ Lq, Strong, Moskalenko, & 
Reimer 2000) is diffuse emission from cosmic ray interac- 
tions with interstellar gas and photons (Hunter et al. 1997). 
However, the LMC is the only external galaxy that has been 
detected in its diffuse y-ray emission (Sreekumar et al. 1992), 
a fact explained by the isotropic flux dilution by distance. At 
1 Mpc, for example, the flux of the Milky Way would approx- 
imately be 2.5 xlO"^ photons cm"^ s"' (>100 MeV), what is 
below the sensitivity achieved up to now in the relevant en- 
ergy domain. The Energetic y-ray Experiment, EGRET, did 
not detect any starburst. Upper limits were imposed for M82, 
F(E > lOOMeV) < 4.4x10"^ photons cm^^ s"', and NGC 253, 
F(E > lOOMeV) < 3.4 x 10"** photons cm"^ s"' (Blom et 
al. 1999), the two nearest starbursts, as well as upon many lu- 
minous infrared galaxies, for which similar constraints were 
found (Torres et al. 2004). These upper limits are barely above 
the theoretical predictions of models for the y-ray emission of 
galaxies, constructed with different levels of detail (see Torres 
2004b for a review). Starbursts and luminous infrared galax- 
ies are expected to compensate the flux dilution produced by 
their relatively larger distance to Earth with their higher cosmic 
ray flux, and become sources for the next generation of instru- 
ments measuring photons in the 100 MeV - 10 GeV regime. 



like the Gamma-ray Large Area Telescope, GLAST (e.g., Volk, 
Aharonian & Breitschwerdt 1996, Paglione et al. 1996, Blom 
et al. 1999, ToiTes et al. 2004, ToiTes 2004). 

In this paper, we analyze one such galaxy, the well stud- 
ied starburst NGC 253. Herein, we present a self-consistent 
multiwavelength modelling of the central region of the galaxy, 
taking into account the latest measurements of densities, su- 
pernova explosion rate, radio emission, and molecular mass, 
among other physical parameters that we use as input for our 
scenario. 

2. CANGAROO observations of NGC 253 

Recently, the CANGAROO collaboration reported the detec- 
tion of NGC 253 at TeV y-ray energies (an llcr confidence 
level claim), observed during a period of two years in 2000 and 
2001 by about 150 hours (Itoh et al. 2002, 2003). Their mea- 
sured differential flux was fitted by Itoh et al. (2003) with a 
power law and an exponential cutoff, obtaining 



— = (2.85 ± 0.71) X 10"'^(£/lTeV) 
dE 



(-3.85.0.46) em-2s-iTeV-\ 
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dF 
"dE 
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with fl = 6 X IQ-^ cm-h-^TeW-\ Eq = 0.0002 TeV, and 
b - 0.25 + 0.01 VTeV. Both parameterizations are sensi- 
ble reproductions of the observational data, although the for- 
mer is clearly prefeiTed for simplicity upon the light of an 
equally good fit. The flux uncertainty are the square root of 
the quadratic sum of the statistical and systematic eiTors. Note 
that the slope of the power law spectrum is very uncertain, but 
steep. Indeed, an extrapolation of this power-law spectrum to 
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lower energies violates the measured upper limits in the GeV 
regime. The CANGAROO collaboration suggested a turn-over 
below the TeV region and proposed the second spectral form. 
They have also claimed that the emission at the highest energies 
is inconsistent with it being produced in a point like source, and 
proposed for it an inverse Compton origin in a kpc-scale y-ray 
halo. 

The HESS array have observed NGC 253 (see below). 
In several other observations of sources that have previ- 
ously been targets for CANGAROO, the HESS collaboration 
have presented results in clear contradiction with the former 
CANGAROO reports. This is most notably the case for SN 
1006 (Aharonian et al. 2005a), PSR B 1706-44 (Aharonian et al. 
2005b), and to some extent also for the supernova remnant RX 
J1713-3857 (Aharonian et al. 2004a), and the Galactic Center 
(Aharonian et al. 2004b). This may suggest some kind of sys- 
tematic difference in the treatment of both sets of observational 
data. Such systematic effect should explain why CANGAROO 
spectra are steeper and their measured fluxes are one order of 
magnitude higher than the upper limits or measurements ob- 
tained with HESS. The CANGAROO collaboration is now cali- 
brating their stereo system, and will be re-observing these prob- 
lematic cases within a year or so (R. Enomoto 2005, private 
communication) . 

In what follows we focus on producing a detailed mul- 
tiwavelength theoretical model for the central region of 
NGC 253, irrespective of CANGAROO measurements (i.e., we 
will not try to reproduce their spectrum, but we will derive pre- 
dictions of fluxes based on a set of well-founded assumptions). 
The aim is to see whether a model based on observations at all 
wavelengths and first principles would -while being consistent 
with multiwavelength testing- predict that the central region of 
NGC 253 alone can produce a sufficiently high y-ray flux so 
as to be detected by the current ground-based Cerenkov tele- 
scopes and future MeV-GeV satellites. The central region of 
the galaxy would look like a point like source for the field of 
view and angular resolution of imaging Cerenkov telescopes. 
Then, we shall explore if we would rather theoretically expect 
a HESS non-confirmation of CANGAROO results regarding 
perhaps both, the flux and the extension. 

3. Phenomenology of the central region of 
NGC 253 

A wealth of new multiwavelength data was obtained for 
NGC 253 during the last decade (after the previous modelling 
by Paghone et al. 1996, which did not include photons energies 
above 100 GeV, see below). It is located at a distance of ~ 2.5 
Mpc (Turner and Ho 1985, Maurbersger et al. 1996) and it is a 
nearly edge-on (incUnation 78°, Pence 1981) barred Sc galaxy. 
The continuum spectrum of NGC 253 peaks in the FIR at about 
100 jum, with a luminosity of 4 x 10'" (Telesco & Harper 
1980, Rice et al. 1988, Melo et al. 2002). The FIR lunoinosity is 
at least a factor of 2 larger than that of our own Galaxy (Cox & 
Mezger 1989, Dudley & Wynn-WiUiams 1999), and it mainly 
(about half of it according to Melo et al.'s (2002) 1 arcmin reso- 
lution ISOPHOT observations) comes from the central nucleus. 



IR emission can be understood as cold {T ~ 50K) dust repro- 
cessing of stellar photon fields. 

When observed at 1 pc resolution, at least 64 individual 
compact radio sources have been detected within the central 
200 pc of the galaxy (Ulvestad & Antonucci 1997), and roughly 
15 of them are within the central arcsec of the strongest ra- 
dio source, considered to be either a buried active nucleus or 
a very compact SNR. Of the strongest 17 sources, about half 
have flat spectra and half have steep spectra. This indicates 
that perhaps half of the individual radio sources are dominated 
by thermal emission from H II regions, and half are optically 
thin synchrotron sources, presumably SNRs. There is no com- 
pelling evidence for any sort of variability in any of the com- 
pact sources over an 8 yr time baseline. The most powerful flat- 
spectrum central radio source is clearly resolved in the study of 
Ulvestad and Antonucci (1997) and appears to be larger than 
the R136 cluster located in 30 Doradus, containing about 10^ 
M0 in stars and 600 M© in ionized gas. The age was estimated 
to be less than 4x10^ yr. The region surrounding the central 
200 pc has also been observed with subarcsec resolution and 22 
additional radio sources stronger than 0.4 mJy were detected 
within 2kpc of the galaxy nucleus (Ulvestad 2000). The region 
outside the central starburst may account for about 20% of the 
star formation of NGC 253, is subject to a supernova explosion 
rate well below 0.1 yr ', and has an average gas density in the 
range 20-200 cm"^, much less than in the most active nuclear 
region of NGC 253 (Ulvestad 2000). 

Carilli (1996) presented low frequency radio continuum ob- 
servations of the nucleus at high spatial resolution. Free-free 
absorption was claimed to be the mechanism producing a flat- 
tening of the synchrotron curve at low energies, with a turnover 
frequency located between 10^'^ and 10^ Hz. The emission 
measures needed for this turnover to happen, for temperatures 
in the order of lO'* K, is at least 10^ pc cm"*". Tingay (2004) 
observed NGC 253 using the Australian Long Baseline Array 
and provided fits with free-free absorption models for the radio 
spectrum of six sources. He concluded that the free-free opac- 
ity in the central region has to be in the range of 1 to 4 at 1.4 
GHz, implying emission measures of a few times 10*" pc cm"*" 
in this particular directions, again for temperatures of the order 
of 10" K. 

As shown by infrared, millimeter, and centimeter observa- 
tions, the 200 pc central region dominates the current star for- 
mation in NGC 253, and is considered as the starburst central 
nucleus (e.g., Ulvestad and Antonucci 1997, Ulvestad 2000). 
Centimeter imaging of this inner starburst, and the limits on 
variability of radio sources, indicates a supernova rate less than 
0.3 yr"' (Ulvestad & Antonucci 1997), which is consistent with 
results ranging from 0. 1 to 0.3 yr~' inferred from models of the 
infrared emission of the entire galaxy (Rieke et al. 1980; Rieke, 
Lebofsky & Walker 1988, Forbes et al. 1993). Van Buren and 
Greenhouse (1994) developed, starting from ChevaUer's (1982) 
model for radio emission from supemovae blast waves expand- 
ing into the ejecta of their precursor stars, a direct relation- 
ship between the FIR luminosity and the rate of supernova ex- 
plosions. The result is "R = 2.3 x 10"'^Lfir/Lo yr"\ which 
is in agreement, within uncertainties, with the previous esti- 
mates. The star formation rate at the central region has been 
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computed from IR observations, resulting in 3.5 Mq yr and 
represents about 70% of the total star formation rate measured 
for NGC 253 (Melo et al. 2002). When compared with Local 
Group Galaxies, the supernova rate in NGC 253 is one order 
of magnitude larger than that of M31, the largest of the Local 
Group (Pavlidou and Fields 2001). 

Paglione et al. (2004) obtained high resolution (5".2 x 5".2) 
interferometric observations of the CO Une 7 = 1 —> in or- 
der to study the structure and kinematics of the molecular gas 
in the central nucleus. This study enhances that of Sorai et al. 
(2000), which, although done with less angular resolution, ob- 
tained compatible results. The general morphology of the CO 
map is consistent with other high resolution studies. It shows 
an extended ridge of emission aUgned with an infrared-bright 
bar and a central group of bright clouds aligned with the ma- 
jor axis of the galaxy, orbiting the radio nucleus in a possible 
ring. The central clouds move around the radio nucleus as a 
solid body, similar to the distribution of the radio sources, cen- 
tral HCN clouds, and central near-infrared emission (Paglione 
et al. 1995, 1997; Ulvestad & Antonucci 1997). Much of the 
molecular gas in NGC 253 appears to be highly excited (Wild 
et al. 1992; Mao et al. 2000; Ward et al. 2003). Observations of 
7 = 4 -> 3 and 7 = 6 -> 5 transitions of CO as well as HCN 
lines suggest the existence of localized spots with values of 
densities in excess of 10"* cm"^ (Israel & Baas 2002, Paglione, 
Jackson, & Ishizuki 1997, PagUone, Tosaki & Jackson 1995, 
Harris et al. 1991). Bradford et al. (2003) have reported CO 
7 = 7^6 observations and also find that the bulk of molecu- 
lar gas in the central 180 pc of the galaxy is highly excited and 
at a temperature of about 120 K. They concluded that the best 
mechanism for heating the gas is cosmic ray bombardment over 
the gas residing in clouds, with density about 4.5 xlO"* cm"^. 

Current estimates of the gas mass in the central 20" - 50" 
(< 600 pc) region range from 2.5 xlO^ M© (Harrison, Henkel 
& Russell 1999) to 4.8 xlO*^ Mq (Houghton et al. 1997), see 
Bradford et al. (2003), Sorai et al. 2000, and Engelbracht et al. 
(1998) for discussions. For example, using the standard CO to 
gas mass conversion, the central molecular mass was estimated 
as 1.8 xlO*^ Mq (Mauersberger et al. 1996). It would be factor 
of ~ 3 lower if such is the correction to the conversion factor 
in starburst regions which are better described as a filled inter- 
cloud medium, as in the case of ULIRGs, instead of a collection 
of separate large molecular clouds, see Solomon et al. (1997), 
Downes & Solomon (1998), and Bryant & Scoville (1999) for 
discussions. Thus we will assume in agreement with the men- 
tioned measurements that within the central 200 pc, a disk of 
70 pc height has ~ 2 xlO^ M© uniformly distributed, with a 
density of ~ 600 cm"^. Additional target gas mass with an av- 
erage density of ~50 cm"^ is assumed to populate the central 
kpc outside the innermost region, but subject to a smaller su- 
pernova explosion rate ~ O.OI yr"\ 10% of that found in the 
most powerful nucleus (Ulvestad 2000). 

The central region of this starburst is packed with massive 
stars. Watson et al. (1996) have discovered four young globular 
clusters near the center of NGC 253; they alone can account 
for a mass well in excess of 1.5xlO®Mo (see also Keto et al. 
1999). Assuming that the star formation rate has been continu- 
ous in the central region for the last lO' yrs, and a Salpeter IMF 



for 0.08-100 Mq, Watson et al. (1996) find that the bolomet- 
ric luminosity of NGC 253 is consistent with 1.5 xlO^M© of 
young stars. Physical, morphological, and kinematic evidence 
for the existence of a galactic superwind has been found for 
NGC 253 (e.g., McCarthy et al. 1987, Heckman et al. 1990, 
Strickland et al. 2000, 2002, Pietsch et al. 2001, Forbes et al. 
2000, Weaver et al. 2002, Sugai, Davies & Ward 2003). This 
superwind creates a cavity of hot (~ 10^ K) gas, with cool- 
ing times longer than the typical expansion timescales. As the 
cavity expands, a strong shock front is formed on the contact 
surface with the cool interstellar medium. Shock interactions 
with low and high density clouds can produce X-ray contin- 
uum and optical line emission, respectively, both of which have 
been directly observed (McCarthy et al. 1987). The shock ve- 
locity can reach thousands of km s ' . This wind has been pro- 
posed as the convector of particles which have been already 
accelerated in individual SNRs, to the larger superwind region, 
where Fermi processes could upgrade their energy up to that 
detected in ultra high energy cosmic rays (Anchordoqui et al. 
1999, Anchordoqui et al. 2003, Torres & Anchordoqui 2004). 



4. Diffuse modelling 

The approach to compute the steady multiwavelength emis- 
sion from NGC 253 follows that implemented in 62-diefuse, 
which we have used with some further improvements (Torres 
2004). The steady state particle distribution is computed within 
(3-DiFFusE as the result of an injection distribution being sub- 
ject to losses and secondary production in the ISM. In gen- 
eral, the injection distribution may be defined to a lesser de- 
gree of uncertainty when compared with the steady state one, 
since the former can be directly Unked to observations. The 
injection proton emissivity, following Bell (1978), is assumed 
to be a power law in proton kinetic energies, with index p, 
2inj(£'p,kin) = ^(£'p,kin/GeV)"^ , where K is a normaUzation 
constant and units are such that [Q]= GeV"' cm"-' s"^ The 
normalization K is obtained from the total power transferred 
by supemovae into CRs kinetic energy within a given volume 

Ie^'ZZ^ 6inj(£'p,kin)£'p,kmC'£'p,km = "-^^p.^n, min/(-P + 2) = 

2, rji'P'Ri/V where it was assumed that p 2, and used the 
fact that /ip, kin mill ^ ^p.kin.max in the second equality. % 
(Zi l^i - is defined as the rate of supernova explosions in the 
star forming region being considered, V being its volume, and 
77,, the transferred fraction of the supernova explosion power 
{r ~ 10^' erg) into CRs (e.g., Torres et al. 2003 and refer- 
ences therein). The summation over i takes into account that 
not all supernovae will transfer the same amount of power into 
CRs (alternatively, that not all supemovae will release the same 
power). The rate of power transfer is assumed to be in the range 
0.05 < rji < 0.15, the average value for 77 being ~ 10%. At low 
energies the distribution of cosmic rays is flatter, e.g., it would 
be given by equation (6) of BeU (1978), correspondingly nor- 
mahzed. We have numerically verified that to neglect this dif- 
ference at low energy does not produce any important change 
in the computation of secondaries, and especially on y-rays at 
the energies of interest. 
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The general diffusion-loss equation is given by (see, e.g., 
Longair 1994, p. 279; Ginzburg & Syrovatskii 1964, p. 296) 

_ D NiE) -^^-^ [bimiE)] - QiE) = -^(1) 
t{E) dE at 

In this equation, D is the scalar diffusion coefficient, Q{E) rep- 
resents the source term appropriate to the production of parti- 
cles with energy E, t{E) stands for the confinement timescale, 
A'(£') is the distribution of particles with energies in the range 
E and E + dE per unit volume, and b(E) - - (dE/dt) is the rate 
of loss of energy. The functions b(E), t(E), and QiE) depend 
on the kind of particles. In the steady state, dN{E)/dt - 0, 
and, under the assumption of a homogeneous distribution of 
sources, the spatial dependence is considered to be irrelevant, 
so that DV^ N(E) = 0. Eq. Q can be solved by using the Green 
function 



G(E,E') = 




such that for any given source function, or emissivity, Q(E), the 
solution is 

dE'Q(E')G(E,E'). (3) 

The confinement timescale will be determined by several con- 
tributions. One on hand, we consider the characteristic escape 
time in the homogeneous diffusion model (Berezinskii et al. 
1990, p. 50-52 and 78) td = R^/(2D(E)) = to/(/3{E /GeWT), 
where /3 is the velocity of the particle in units of c, R is the 
spatial extent of the region from where particles diffuse away, 
and D(E) is the energy-dependent diffusion coefficient, whose 
dependence is assumed oc E*^, with jj. ~ 0.5. ' tq is the char- 
acteristic diffusive escape time at ~ 1 GeV. tq for NGC253 is 
not known. One can only assume its value and compare it with 
that for other galaxies (e.g. our own Galaxy, or M33, Duric et 
al. 1995); the value we choose also parallels that obtained in 
an earUer study of NGC 253 or on M82 (PagHone et al. 1996, 
Blom et al. 1999). We analyze the sensitivity of the model to 
To below. On the other, the total escape timescale will also take 
into account that particles can be carried away by the collec- 
tive effect of stellar winds and supernovae. t^, the convective 
timescale, is ~ R/V, where V is the collective wind veloc- 
ity. For a wind velocity of 300 km/s and a radius of about 
the size of the innermost starburst (see below), this timescale 
is less than a million years ( 3x10^ yr). The outflow veloc- 
ity is not very well known, however, but minimum reasonable 
values between 300 and 600 km s ' have been claimed, and 
could even reach values of the order of thousand of km s ' 
(Strickland et al. 2002). Pion losses (which are catastrophic, 
since the inelasticity of the collision is about 50%) produce 
a loss timescale = {dEjdty™ jE (see, e.g., Mannheim & 

' We emphasize that the use of an homogeneous model is an as- 
sumption, but justified in the compactness of the innermost starburst 
region. We are basically assuming that there is an homogeneous distri- 
bution of supernovae in the central hundreds of pc, what is supported 
observationally (Ulvestad and Antonucci 1997). 



Schlickeiser 1994), which is similar in magnitude to the con- 
vective timescale and dominates with it the shaping of the pro- 
ton spectrum. Thus, in general, for energies higher than the 
pion production threshold t"' (£) = t^' -H t^"' -H Tpp' . For elec- 
trons, the total rate of energy loss considered is given by the 
sum of that involving ionization, inverse Compton scattering, 
bremsstrahlung, and synchrotron radiation. We have also in- 
corporated adiabatic losses. Full Klein-Nishina cross section is 
used while computing photon emission, and either Thomson 
or extreme Klein-Nishina approximations, as needed, are used 
while computing losses. For the production of secondary elec- 
trons, Q-DiFFUSE computes knock-on electrons and charged pion 
processes producing both electrons and positrons. In the case of 
y-ray photons, we compute bremsstrahlung, inverse Compton 
and neutral pion decay processes. For the latter, an Appendix 
provides a more detailed discussion of the different approaches 
to compute the neutral pion-induced y-ray emissivity. For radio 
photons, we compute, using the steady distribution of electrons, 
the synchrotron, and free-free emission. Free-free absorption is 
also considered in order to reproduce the radio data at low fre- 
quencies. The FIR emission is modelled with a dust emissivity 
law given by v^Bie, T), where cr = 1.5 and B is the Planck 
function. The computed FIR photon density is used as a target 
for inverse Compton process as well as to give account of losses 
in the y-ray scape. The latter basically comes from the opacity 
to yy pair production with the photon field of the galaxy nu- 
cleus. The fact that the dust within the starburst reprocesses the 
UV star radiation to the less energetic infrared photons implies 
that the opacity to yy process is significant only at the highest 
energies. The opacity to pair production from the interaction of 
a y-ray photon in the presence of a nucleus of charge Z is con- 
sidered too. For further details and relevant formulae see Torres 
(2004). 

4.1. Comparison witii previous models 

When compared with the previous study on high energy 
emission from NGC 253, by Paglione et al. (1996), several 
methodological and modelling differences are to be mentioned. 
Paglione et al.'s assumed distance, size, gas mass, density, and 
supernova explosion rate of the central region are different from 
those quoted in Table 1 . The former authors modelled, based on 
earlier data (e.g., Canzian et al. 1988) a starburst region at 3.4 
Mpc (a factor of 1.36 farther than the one currently adopted), 
of 325 pc radius (about 3 times larger than the one adopted 
here). This region is larger than what is implied by the cur- 
rent knowledge of the central starburst, where the supernova 
explosion rate Paglione et al. used is actually found and the 
cosmic ray density is maximally enhanced. The average den- 
sity assumed by Paglione et al., 300 cm"-', gives a target mass 
~ 2 X 10*^ Mq, which is at the upper end of all current claims 
for the central nucleus, or already found as excessive. The tar- 
get mass of the innermost region differs from ours by a factor 
of about 6, ours being smaller. The fraction of the supernova 
explosion converted into cosmic rays (20% for Paglione et al., 
a factor of 2 larger than ours) seems also excessive in regards 
of the current measurements of SNR at the highest energies. 
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We have also considered, especially to test its influence in the 
total y-ray output, a surrounding disk with a smaller supernova 
rate, following the discovery of several SNRs in that region 
(Ulvestad 2000). Finally, Paglione et al. (1996) study did not 
produce results above 200 GeV. ^ 

The (3-DiFFUSE set uses different parameterizations for pion 
cross sections as compared with those used by Marscher and 
Brown (1978), whose code was the basis of Paglione et al.'s 
study. Our computation of neutral pion decay y-rays is dis- 
cussed in the Appendix. In addition, (3-diffuse uses the full 
inverse Compton Klein-Nishina cross section, computes sec- 
ondaries (e.g., knock-on electrons) without resorting to param- 
eterizations which are valid only for Earth-like cosmic ray (CR) 
intensities, fixes the photon target for Compton scattering start- 
ing from modelling of the observations in the FIR, and consid- 
ers opacities to y-ray scape. 

5. Results 

5.1. Summary of model parameters 

We begin the discussion of our results by presenting a sum- 
mary of the parameters we have used for, and obtained from, 
our modelling. These are given in Table 1 . There, the mark OM 
refers to Obtained from modelling and ST or see text refers 
to parameters discussed in more detail in the previous Section 
on phenomenology, where references are also given. These pa- 
rameters values or ranges of values are fixed by observations. 
Finally, the mark A refers to assumed parameters, in general 
within a range. Variations to the values given in Table 1 are 
discussed below. 

5.2. Steady proton and electron population 

The numerical solution of the diffiision-loss equation for pro- 
tons and electrons, each subject to the losses previously de- 
scribed, is shown in FigureQ] We have adopted a diffusive resi- 
dence timescale of 10 Myr, a convective timescale of 1 Myr, 
and a density of ~ 600 cm"^. In the case of electrons, the 
magnetic field with which synchrotron losses are computed in 
FigureQ]is 300 fiG. The latter is fixed below, requiring that the 
steady electron population produces a flux level of radio emis- 
sion matching observations. An injection electron spectrum is 
considered -in addition to the secondaries- in generating the 
steady electron distribution. The primary electron spectrum is 
assumed as that of the protons times a scaling factor; the in- 
verse of the ratio between the number of protons and electrons, 
Np/Ne (e.g., BeU 1978). This ratio is about 100 for the Galaxy, 
but could be smaller in star forming regions, where there are 
multiple acceleration sites. For instance, Volk et al. (1989) ob- 

^ To further ease the comparison, we here note some typos in 
Paglione et al. (1996) paper: The factor b(E) should be elevated to the 
minus one in their equation (4), as well as the term in their equation 
(3). The y-axis of their Figure 1 is not the emissivity, but the emis- 
sivity divided by the density; units need to be accordingly modified, 
see e.g. Abraham et al. (1966). Ep in their equation (7) and x-axis of 
Figure 2 and 3 is the kinetic energy, but the generic E in Equation (1) 
is the total energy. The y-axis of Figure 2 is in units of cm"^ GeV"' . 



tained Np/Ne ~ 30 for M82. Np/Ne - 50 is assumed for the 
central disk of NGC 253. From about Ee - ~ 10"' to ~ 10 
GeV, the secondary population of electrons (slightly) domi- 
nates, in any case. This is shown in Figure |5] It is in this re- 
gion of energies where most of the synchrotron radio emission 
is generated, and thus the ability of producing a high energy 
model compatible with radio observations is a cross check for 
the primary proton distribution. 

Figure|2lshows the ratio of the steady proton population in 
the SD to that in the IS. Because the SD is subject to a smaller 
supernova explosion rate, it has an smaller number of protons 
in its steady distribution, at all energies, of the order of 1% 
of that of the IS. Then, it will play a subdominant role in the 
generation of y-ray emission, as we show below. In the right 
panel of Figure |2 and for further discussion in the following 
Sections, we present the ratio between the steady proton distri- 
bution in the IS, when the gas density is artificially enhanced 
and diminished by a factor of 2 from the assumed value of 600 
cm"-'. 

5.3. IR and radio emission 

The continuum emission from NGC 253, at wavelengths be- 
tween ~ 1 cm and ~ 10 microns, was measured by several au- 
thors, e.g., see Figurel^and Section 3. These observations did 
not in general distinguish, due to angular resolution, only the 
emission coming from the innermost starburst region. Instead, 
they also contain a contribution coming from photons produced 
in the surrounding disk and farther away from the nucleus. 
The IR continuum emission is mainly produced thermally, by 
dust, and thus it could be modelled with a spectrum having 
a dilute blackbody (graybody) emissivity law, proportional to 
v^B{e, T), where B is the Planck function. Figurel^shows the 
result of this modelling and its agreement with observational 
data when the dust emissivity index cr = 1 .5 and the dust tem- 
perature Tdust = 5QK. Diff'erent total luminosities, the normal- 
ization of the dust emission (see the appendix of Torres 2004 
for details), are shown in the Figure to give an idea of the con- 
tribution of the innermost region with respect to that of the rest 
of the galaxy. According to Melo et al. (2002), about half of 
the total IR luminosity is produced in the IS, what is consis- 
tent with the data points being intermediate between the curves 
with LiR 2x10'° and 4x10'° Lq, since the latter were obtained 
with beamsizes of about 20-50 arcsec (~ 240-600 pc at the 
NGC 253 distance). 

The influence of the magnetic field upon the steady state 
electron distribution is twofold. On one hand, the greater the 
field, the larger the synchrotron losses -what is particularly 
visible at high energies, where synchrotron losses play a rel- 
evant role. On the other, the larger the field the smaller the 
steady distribution. These effects evidently compete between 
each other in determining the final radio ffux. The magnetic 
field is required to be such that the radio emission generated by 
the steady electron distribution is in agreement with the obser- 
vational radio data. This is achieved by iterating the feedback 
between the choice of magnetic field, the determination of the 
steady distribution, and the computation of radio ffux, addition- 
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Table 1. Measured, assumed, and derived values for different physical quantities at the innermost starburst region of NGC 253 
(IS), a cylindrical disk with height 70 pc, and its surrounding disk (SD). 



Physical parEinctcrs 


Svmbol 


Value 


Units 


Comment 


111 cttiTir*!^ 


J) 




Mpc 




Tnr*l in citir^n 
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ST 
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2 X 10'" 




ST 


Radius of the IS 




100 




ST 


Radius surrounding disk (SD) 




1000 


DC 


ST 


TTnifViTTn HpTTsitv of the T.S 


"IS 
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-3 

cm 


ST 
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ST 
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Mr. 


ST 


Gas mass of the SD 


-'"SD 
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Mr-. 


ST 


Qiinprtinvn pvnloQion ratp of trip 

ij LtlJtl IIU V£l v-AUHJolVJll IClLl^ \Jl Lilt- AvJ 


^ 


~ 0.08 


SN yr"' 


ST 


Supernova explosion rate of the SD 




~ 0.01 


SN yr"' 


ST 


Tvnipal *;iinpt*rinvji PTmlnQinn ptiptov 




1051 


erg 


ST 


SN energy transferred to cosmic rays 


77 


~ 10 


% 


ST 
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y 


300 - 600 


km s~' 


ST 
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(T 


1 s 




OM 


r^iiist tpmnpTiitiiT'p 


T 

^ dust 


50 


K 


OM 


Emission measure 


EM 
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pc cm~^ 


OM 


Ionized gas temperature 


T 


10^ 


K 


OM 


Magnetic field of the IS 


B 


300 


juG 


OM 


Slope of primary injection spectrum 


P 


2.2-2.3 




A 


Maximum energy considered for primaries 




100 


TeV 


A 


Diffusion coefficient slope 


ft 


0.5 




A 


Proton to electron primary ratio 




50 




A 


Diffusive timescale 


To 


1-10 


Myr 


A 




Fig. 1. Steady proton (left panel) and electron (right panel) distributions in the innermost region of NGC 253. 



ally taking into account free- free emission and absorption pro- 
cesses. Whereas free-free emission is subdominant when com- 
pared with the synchrotron flux density, free-free absorption 
plays a key role at low frequencies, determining the opacity. 
We have found a reasonable agreement with all observational 
data for a magnetic field in the innermost region of 300 jiG, an 
ionized gas temperature of about 10"* K, and an emission mea- 
sure of 5 X 10^ pc cm~^, the latter two are in separate agreement 



with the free-free modelhng of the opacity of particular radio 
sources, as discussed in Section 3. The value of magnetic field 
we have found for the IS is very similar to that found for the 
disk of Arp 220 (Torres 2004) and compatible with measure- 
ments in molecular clouds (Crutcher 1988, 1994, 1999). 
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Fig. 2. Left: Ratio of the steady proton population in the surrounding disk to that in the innermost starburst region. Right: Ratio 
between the steady proton distribution in the IS, when the gas density is artificially enhanced and diminished by a factor of 2. 
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Fig. 3. Steady population of primary-only and secondary-only electrons. Only the region of the secondary dominance of the 
distribution is shown. 



5.4. y-ray emission 



In the left panel of Figure|5] bremsstrahlung, inverse Compton, 
and pion decay y-ray fluxes of the IS are shown together with 
the total contribution of the SD and the total diff'erential flux 
of the whole system. These results are obtained with the model 
just shown to be in agreement with radio and IR observations. 
As mentioned before, the contribution of the SD, even when 
having a factor of ~ 8 more mass than the IS, is subdominant. 
The reason for this needs to be found in that this region is much 
less active (Ulvestad 2000). 



Our predictions, while complying with EGRET upper lim- 
its, are barely below them. If this model is correct, NGC 253 is 
bound to be a bright y-ray source for GLAST. 

The integral fluxes are shown in the right panel of Figure|5] 
Our model complies again with the integral EGRET upper limit 
for photons above 100 MeV, and predicts that, given enough 
observation time, NGC 253 is also to appear as a point-like 
source in an instrument like HESS. Note, however, that quite 
a long exposure may be needed to detect the galaxy, and also, 
that our fluxes are only a few percent of those reported by the 
CANGAROO coUaboration. 
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Fig. 4. Left: IR flux from NGC 253 assuming a dilute blackbody with temperature Jdust = 50 K and diff"erent total luminosities. 
Right: Multifrequency spectrum of NGC 253 from radio to IR, with the result of our modelling. The experimental data points 

correspond to: pentagons, Melo et al. (2002); diamonds, Telesco et al. (1980); down-facing triangles, Rieke et al. (1973); stars, 
Hildebrand et al. (1977); up-facing triangles, Ehas et al. (1978); circles, Ott et al. (2005); squares, Carilh (1996). 




Fig. 5. Left: Differential y-ray fluxes from the central region of NGC 253. Total contribution of the surrounding disk is separately 
shown, as are the EGRET upper limits. In the case of the IS, we separately show the relative contributions of bremsstrahlung, 
inverse Compton, and neutral pion decay to the y-ray flux. Right: Integral y-ray fluxes. The EGRET upper limit (for energies 
above 100 MeV), the CANGAROO integral flux as estimated from their fit, and the HESS sensitivity (for a 5cr detection in 50 
hours) are given. Absorption effects are already taken into account. Also shown is the recently released HESS upper limit curve 
on NGC 253. 



An additional source of TeV photons not considered here is 
the hadronic production in the winds of massive stars (Romero 
& Torres 2003). However, a strong star forming region such 
as the nucleus of NGC 253 is bound to possess much more 
free gas than that contained within the winds of massive stars, 



which albeit numerous, have mass loss rates typically in the 
range of IQ-^-lQ-' M^? 

^ In Romero & Torres (2003), higher mass loss rates up to 10"' Mg, 
i.e., grammages between 50 and 150 g cm"^ were allowed. These val- 
ues, although have been found in perhaps one or two Galactic early O 
stars, are uncommon. Since the size of the base of the wind for each 
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Fig. 6. Left: Opacities to y-ray scape as a function of energy. The highest energy is dominated by yy processes, whereas yZ 
dominates the opacity at low energies. Significant Tmax are only encountered above 1 TeV. Right: Modification of the y-ray 
spectrum introduced by the opacity to y-ray escape. 



5.5. Opacities to y-ray escape 

Two sources of opacities are considered: pair production from 
the y-rays interaction with the photon field or with the charged 
nucleus present in the medium. Compton scattering and atten- 
uation in the magnetic field by one-photon pair production are 
negligible. 

The opacity to yy pair production with the pho- 
ton field, which, at the same time, is target for 
inverse Compton processes, can be computed as 

poo 

(/i/cos(/)) n(e)cre-e+ie,Ey)de where e is the energy 
of the target photons, Ey is the energy of the y-ray 
in consideration, Rc is the place where the y-ray pho- 
ton was created within the system, and o"^ ,;+(e, Zs^)^'' = 
(3crW16)(l -yS2)(2yS(^2 _ 2) + (3 _^)))^ ^jth 

/3 - (1 - (mc^)^ /(e Ey)y^^ and ctj- being the Thomson cross 
section, is the cross section for yy pair production (e.g. Cox 
1999, p. 214). Note that the lower limit of the integral on e in 
the expression for the opacity is determined from the condition 
that the center of mass energy of the two colliding photons 
should be such that /3 > 0. The fact that the dust within the 
starburst reprocesses the UV star radiation to the less energetic 
infrared photons implies that the opacities to yy process is 
significant only at the highest energies. No source of this kind 
of opacity is assumed outside the system under consideration, 
since the nearness of NGC 253 makes negligible the opacity 
generated by cosmological fields. 



star, the grammage, and the ambient enhancement of cosmic rays were 
independently allowed to take values within their assumed ranges in 
the Monte Carlo simulation of Romero & Torres (2003), the stars with 
the most favorable parameters would dominate the sum, overestimat- 
ing the relative importance of their fluxes. 



The cross section for pair production from the interac- 
tion of a y-ray photon in the presence of a nucleus of charge 
Z in the completely screened regime (Ey/mc^ » l/{aZ)) 
is independent of energy, and is given by (e.g. Cox 1999, 
p.213) cr% = (3ci'ZVr/2;T)(7/91n(183/Zi/3) _ i/^^y 

lower energies the relevant cross section is that of the no- 
screening case, which has a logarithmic dependency on energy, 
crl^^^ = (3afZVr/27r)(7/91n(2£y/mc2)- 109/54), and matches 
the complete screening cross section at around 0.5 GeV. Both 
of these expression are used to compute the opacity, depending 
on Ey. 

In the left panel of Figure |5] we show the different contri- 
butions to the opacity. The equation of radiation transport ap- 
propriate for a disk is used to compute the predicted y-ray flux 
taking into account all absorption processes (see Appendix of 
Torres 2004 for details). The right panel of Figure|6lshows the 
eff'ect of the opacity on the integral y-ray fluxes, only evident 
above 3 TeV. 

5.6. Exploring tine parameter space and degeneracies 

As it is summarized in Table 1, most of the model parameters 
are well fixed from observations. There are however a couple 
of assumptions which, whereas being not well bounded from 
experiments, may produce slight degeneracies. Consider for in- 
stance the proton injection slope p and the diff'usive scale tq. 
For the former we have assumed p - 2.3, which is in agree- 
ment with the recent results from HESS regarding y-ray ob- 
servations at TeV energies of supernova remnants and uniden- 
tified extended sources. However, it would be certainly within 
what one would expect from proton acceleration in supernova 
remnants, and also within experimental uncertainty, if a better 
description for the average proton injection slope in NGC 253 
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is 2.2 instead of 2.3. Table |2 shows the influence of this kind 
of choice on our final results. A harder slope slightly increases 
the integral flux. Similarly, the diffusive timescale is not well 
determined, and it may be arguable perhaps within one order 
of magnitude. Table|2also shows the influence of this choice. 
Ultimately, high energy y-ray observations (from GeV to TeV) 
are the ones to impose constraints on these values. In any case, 
we remark that pp interaction and convection timescales are 
much shorter (< 1 Myr) and thus dominate the form of N{E). 
To show this in greater detail we show in Figure the result 
for the proton distribution when different convective and the pp 
timescales are taken into account as compared with the solution 
when t{E) - tq, i.e., diffusion only. Clearly, convection plus 
pp timescales dominates the spectrum. 

Regarding degeneracies, both the proton slope and the con- 
finement timescales, however, cannot be much different from 
what we have assumed. If the former were to differ signif- 
icantly, it would be impossible to reproduce the radio data, 
which is the result of the synchrotron emission of the secondary 
electrons. Changes in the number of protons in the IS would 
imply a change in the magnetic field to reproduce radio obser- 
vations, what clearly cannot be pushed much either 

In a less impacting way, varying the value of Np/Ng can 
also vary the results. This variation would be slight because 
of the influence of the more numerous secondary electrons in 
the energetic region of interest for radio emission. On the same 
track, varying the diffusion coefficient fi does not import sub- 
stantial changes. And if the maximum proton energy were to 
differ from the value of 100 TeV we have assumed (what we 
do not expect it to happen in a significant way, since we do 
now observationally know that supernova remnants are sources 
of ~ 10 TeV photons), the end of the spectrum would slightly 
shift accordingly. 

Even within an artificially enlarged uncertainty of the gas 
density, the results will not be modified much: if for any reason 
the average particle density were to be a factor of 2 smaller or 
larger, the y-ray integral flux variations would be within 4% for 
energies above 100 MeV, and within 25% for energies above 
200 GeV. Table|3]shows these results by presenting the integral 
fluxes above a given threshold if the assumed density of 600 
cm"^ is doubled or halved. As can be seen in the right panel 
of Figure |2l if the density is larger (smaller) by a factor ~ 2, 
the resultant steady proton distribution from the same proton 
injection population is smaller (bigger) by a similar factor over 
a wide range of proton energies. As y-ray emissivities are pro- 
portional to both the medium density and the number of steady 
protons, the variations in y-ray fluxes are greatly compensated. 

5. 7. Energetics and cosmic ray enliancement 

The left panel of Figure |8] presents the energy density con- 
tained in the steady proton population above a certain en- 
ergy, i.e., based on Figure [2 the curve shows the integral 
J^Np(Ep)EpdEp. The total energy density contained by the 
steady population of cosmic rays above 1 GeV is about lO"-' of 
the power emitted by all supernova explosions in the last 5 mil- 



lion years. The energy density contained in the steady electron 
population is orders of magnitude less important. 

The cosmic ray enhancement is a useful parameter in es- 
timations of y-ray luminosities in different scenarios. It is de- 
fined as the increase in the cosmic ray energy density with re- 
spect to the local value, cocr,0(E) = Nps{Ep)EpdEp, where 
Npd, is the local cosmic ray distribution obtained from the mea- 
sured cosmic ray flux that we quote in the Appendix. The 
enhancement factor g is then a function of energy given by 
g(E) - ( J^Np(Ep) EpdEp)/a)cR,Q(E). Values of enhancement 
for NGC 253 were proposed g < 3000 for energies above 1 
GeV (e.g., Suchkov et al. 1993), and we can actually verify 
this in our model. The right panel of Figure|S]presents the en- 
hancement factor as a function of proton energy. The larger the 
energy, the larger the enhancement, due to the steep decline 
(oc E^^-''^) of the local cosmic ray spectrum. 

6. HESS observations 

The HESS array has just released (Aharonian et al. 2005c) their 
results for NGC 253. These are based on data taken during the 
construction of the array with 2 and 3 telescopes operating. The 
total observation time was 28 hs, with a mean zenith angle of 
about 14 degrees. Only events where at least two telescopes 
were triggered were used, to enable stereoscopic reconstruc- 
tion. The energy threshold for this dataset was 190 GeV. Upper 
limits from H.E.S.S. on the integral flux of y-rays from NGC 
253 (99 % confidence level) are shown, together with our pre- 
dictions, in the right panel of Figure |5] and zoomed in the re- 
gion above 100 GeV in Figure |5] As an example, above 300 
GeV, the upper limit is 1.9 x 10 photons cm"^ s It can 
be seen that our predictions are below these upper limits at all 
energies but still above HESS sensitivity for reasonable obser- 
vation times. 

7. Concluding remarks 

We have presented a multifrequency model of the central re- 
gion of NGC 253. Following recent observations, we have 
modelled an innermost starburst with a radius of 100 pc and 
a supernova explosion rate of 0.08 yr"', and a surrounding disk 
up to a 1 kpc in radius with an explosion rate about tenfold 
smaller. As a result of our modelling we have found that a 
magnetic field of 300 fiG for the innermost region is consistent 
with high resolution radio observations, with the radiation at 
1 GHz being mostly produced by secondary electrons of cos- 
mic ray interactions. The magnetic field found for the inner- 
most part of NGC 253 is typical of dense molecular clouds 
in our Galaxy, and is close to the (270 fiG) value proposed 
by Weaver et al. (2002) using the equipartition argument. We 
have estimated free-free emission and absorption, and consid- 
ered opacities to the y-ray escape. The hard X-ray emission 
from IC and bremsstrahlung processes produced in this model 
is below observational constraints, e.g., by OSSE, in agreement 
with previous estimations of bremsstrahlung diffuse emission 
(Bhattacharya et al. 1994). This is consistent with measure- 
ments in the center of the Galaxy, where INTEGRAL have 
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Fig. 7. Proton distribution when different convective and the pp timescales are taken into account as compared with a (arbitrary) 
solution when riE) = td, i.e., diffusion only. Clearly, convection plus pp timescales dominates the spectrum. 

Table 2. Exploring the parameter space for p and tq. The results of our adopted model are given in the first column. These results 
already take into account the opacity to photon escape. 



F(E > Eo) 


p = 23 


p = 2.3 


p = 2.3 


p = 2.2 


p = 2.2 


p = 2.2 




To = 10 Myr 


To = 5 Myr 


To = 1 Myr 


To = 10 Myr 


To = 5 Myr 


To = 1 Myr 


Eo = 100 MeV 


2.32E-8 


2.36E-8 


2.21E-8 


2.95E-8 


2.97E-8 


2.75E-8 


Eo = 200 GeV 


1.60E-12 


1.23E-12 


4.76E-13 


4.04E-12 


3.08E-12 


1.15E-12 


Eo = 600 GeV 


3.61E-13 


2.67E-13 


8.98E-14 


l.OOE-12 


7.34E-13 


2.40E-13 


Eo = 1 TeV 


1.78E-13 


1.29E-13 


4.10E-14 


5.16E-13 


3.70E-13 


1.14E-13 


Eo = 2 TeV 


6.29E-14 


4.46E-14 


1.31E-14 


1.92E-13 


1.35E-13 


3.87E-14 



Table 3. The effect of the medium gas density on the y-ray integral fluxes. Results provided are in units of photons cm ^ s \ and 
already take into account the opacity to photon escape. 



F(E > Eo) 


n=300 cm"^ 


«=600 cm" ' 


n=1200cm-^ 


£o = 100 Mev 


2.22E-8 


2.32E-8 


2.37E-8 


Eo = 200 Gev 


1.20E-12 


1.60E-12 


1.95E-12 


Eo = 600 Gev 


2.60E-13 


3.61E-13 


4.52E-13 


Eo = 1 Tev 


1.26E-13 


1.78E-13 


2.26E-13 


Eg = 2 Tev 


4.36E-14 


6.29E-14 


8.09E-14 



shown that hard X-ray emission is not diffusively, but produced 
by point like sources (Lebrun et al. 2004). 

The flux predicted is based on a set of a few well founded 
assumptions, mainly a) that supernova remnants accelerate 
most of the cosmic rays in the central region of NGC 253, and 
b) that they interact with the present gas, whose amount has 
been measured using a variety of techniques. The low opac- 
ity to y-ray escape secure that basically all y-rays produced in 
the direction towards Earth reach us. Observational constraints 
establishes the values of the supernova explosion rate and gas 
content (see Section 2 for references). 

The ease of all the assumptions made in our model, its con- 
currence with all observational constraints, and the unavoid- 



ability of the processes analyzed, lead us to conclude that 
1) GLAST win detect NGC 253, being our predicted lumi- 
nosity (2.3 X 10"^ photons cm"^ s ' above 100 MeV) weU 
above its 1 yr all sky survey sensitivity (GLAST Science 
Requirements Doc. 2003); 2) that our predicted TeV fluxes are 
about one order of magnitude smaller than what was claimed 
by CANGAROO, and thus, that perhaps in this case, a similar 
problem to that found in other sources affected their data tak- 
ing or analysis, and 3) that HESS could detect the galaxy as a 
point like source provided it is observed long enough with the 
full array (> 50 hours, for a detection between 300 and 1000 
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Fig. 8. Left: Energy density contained in tiie steady population of protons above the energy given by the x-axis. Right: Cosmic 
ray enhancement factor obtained from the steady spectrum distribution in the innermost starburst nucleus of NGC 253. 
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Fig. 9. Integral y-ray fluxes zoomed in the region above 100 GeV, together with the (2 telescopes) HESS upper limits and the (4 
telescopes) HESS sensitivity. 



GeV.) We finally note that this model predicts a steady y-ray 
source, so that a posteriori variability estimators (e.g., Torres et 
al. 2001) can be checked for consistency. 



HESS site latitude provides that NGC 253 can be observed very 
close to the zenith (the minimum zenith angle for NGC 253 from 
HESS site is 2 degrees). As a consequence, HESS observations of 
NGC 253 can be done with the minimum energy tlireshold of the ex- 
periment. The MAGIC Telescope, although being at a northern hemi- 
sphere site, is also able to observe NGC 253 at a larger zenith angle, 
about 53 degrees. 



Appendix: Parameterizations of proton-proton 
cross sections for neutral pion decay 

The emissivity resulting from an isotropic intensity of pro- 
tons, Jp(Ep), interacting with fixed target nuclei with number 



density n, through the reaction p+p 
is given by (e.g., Stecker 1971) 



p+p+n — > p+p+2y. 



f 

Eth 



dEp Jp(Ep) 



d(T(E„i> , Ep) 
dEjM 



(4) 



where E^"" 



is the maximum energy of protons in the system, 
and E,h{E^) is the minimum proton energy required to pro- 
duce a pion with total energy E^, and is determined through 
kinematical considerations. It is obtained using the invariant, 

^fs = {liUpiEp + nipij'^ = {M\ + E*J--miyi^ + E*„^, where jis 
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the square of the total energy in the center-of-mass system, Mx 
depends on the reaction channel and represents the invariant 
mass of the system consisting of all particles except the pion, 
E* is the CMS energy of the produced meson (that is connected 
with the laboratory system energy via a Lorentz transformation, 
see Appendices of Moskalenko & Strong 1998 and Blattnig et 
al. 2000b). = (i - + ml)/(2 ^/s), so that for a given value 
of 5, E* will be maximum when takes its minimum value. 
For the case of neutral pion production - Inip. The labo- 
ratory system pion energy, obtained from £*, can be put as a 
function of s. Inverting this relation, thus obtaining s - s{E„), 
and use of s - s{E„) - 2mp(Ep + nip) allow for the minimum 
proton energy to be derived. Finally, do-{E„o,Ep)/dE„o is the 
differential cross section for the production of a pion with en- 
ergy E^o, in the lab frame, due to a collision of a proton of 
energy Ep with a hydrogen atom at rest. 

The y-ray emissivity is obtained from the neutral pion 
emissivity Q„o as 



QyiEy). 



dE. 



EJiE,) 



" (Ei,-miyyn 



(5) 



where E'™\Ey) - Ey + rrr„c /(4Ey) is the minimum pion en- 
ergy required to produce a photon of energy Ey (e.g., Stecker 
1971), and E'^"^(E"j"^) is the maximum pion energy that the 
population of protons can produce. It is obtained using the in- 
variant equation for the maximum proton energy resident in the 
system. 

Thus, an accurate knowledge of the differential cross sec- 
tion for pion production becomes very important to estimate the 
y-ray emissivity. Note that dcr{E„o,Ep)/dE„o can be thought 
as containing the inclusive total inelastic cross section (i.e. 
the cross section multiplied by the average pion multiplicity). 
This can be stated explicitly as done, for instance, in Dermer's 
(1986a) equation 3. 

6-function approximation 

In this formalism (Aharonian and Atoyan 2000), 

Q„o{E„o) = 47m I dEp Jp{Ep) 6{E„o - KE^^n) cr{Ep) , 

Ann I 2 Ejfi \ I , ^;r» \ 

^—Jp\mpC + —\a-\mpC + — I (6) 

where cr is the total cross-section of inelastic pp collisions, and 
K is the mean fraction of the kinetic energy fikin - Ep- nipC^ of 
the proton transferred to the secondary meson per collision. In 
a broad region from GeV to TeV energies, k ~ 0.17. In this ap- 
proximation, then, an accurate knowledge of the total inelastic 
cross section is needed to compute the y-ray emissivity. 

Aharonian and Atoyan (2000) proposed that, since from the 
threshold at isan ~ 0.3 GeV the cross section appears to rise 
rapidly to about 30 mb at energies about Ei^„ ~ 2 GeV, and 
since after that energy it increases only logarithmically, a suffi- 
ciently good approximation is to assume 

o- ~ 30 (0.95 + 0.061n(£kin/GeV))mb for E > IGeV 
~ otherwise. (7) 



It can be seen (e.g., Dermer 1986a) that different parameteriza- 
tions of the cross section below 1 GeV do not noticeably affect 
the results of y-ray emissitivities, since most of the y-rays are 
generated by primary protons having more energy than a few 
GeV, provided the spectrum of primaries is sufficiently broad. 

In a recent paper, Kamae et al. (2005) introduced the effect 
of diffractive interactions and scaling violations in pp — > tt" 
interactions. The diffractive interactions contribution was usu- 
ally neglected in all other computations of y-ray emissivity 
from neutral pion decay to date, and thus one would expect 
an increase in the predicted fluxes. Kamae et al.'s best model, 
dubbed A, for the inelastic (not inclusive) cross section is given 
in Table 1 of their paper, columns 2 and 3. When one compares 
the sum of both diffractive and non-diffractive contributions of 
Kamae et al.'s model with the Aharonian and Atoyan's for- 
mula, one sees that the latter produces an actually larger (but 
quite close) cross section. FigureQ^lshows these results above 
proton kinetic energies of 1 GeV, as well as the difference be- 
tween these cross sections. Kamae et al.'s model A was com- 
pared with Hagiwara's (2002) compilation of pp cross section 
measurements and found in good agreement. When multiplic- 
ity is taken into account, Kamae et al.'s model also agrees with 
the data on inclusive cross sections, a point we discuss in more 
detail below. 

Figure ^2 shows a comparison of the y-ray emissivity ob- 
tained when using Kamae et al.'s model A and equation (|6}. 
Curves are practically indistinguishable in this scale, and their 
ratio is well within a factor of ~ 1.3. In this comparison, the 
proton spectrum is the Earth-like one, Jp{Ep) - 2.2 E'^^-''^ pro- 
tons cm"^ s"' sr"' GeV"' and n - I cnT^. The resulting y-ray 
emissivity is multiplied by 1 .45 to give account of the contribu- 
tion to the pion spectrum produced in interactions with heav- 
ier nuclei both as targets and as projectiles (Dermer 1986a). 
Aharonian and Atoyan's expression for the cross section pro- 
duces a slightly larger value of emissivity than that obtained 
with Kamae's model A, including non-diffractive interactions. 

Differential cross sections parameterizations 

Recently, Blattnig et al. (2000) developed parameterizations of 
the differential cross sections. They have presented a param- 
eterization of the Stephens and Badhwar's (1981) model by 
numerically integrating the Lorentz-invariant differential cross 
section (LIDCS). The expression of such parameterization is 
divided into two regions, depending on the (laboratory frame) 
proton energy (Blattnig et al. 2000, see their equations 23 and 
24). Blattnig et al. have also developed an alternative parame- 
terization that has a much simpler analytical form. It is given 
by 



do-(E„i>,Ep) 

dEjfi 
with 

A = ( -5.8 - 



= e'^mb GeV"' 
1.82 



(8) 



13.5 



4.5 



{Ep - nipf -^ (E„o - »v)"-2 (£,0 - »v)0-4 / 

Both parameterizations were integrated and compared with ex- 
perimental results up to ~ 50 GeV in Blattnig et al.'s (2000) 
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Fig. 10. Comparison between Kamae et al.'s model A, sum of difFractive and non-difFractive contributions with Aharonian and 
Atoyan's formula for the total inelastic cross section. 
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Fig. 11. Comparing y-ray emissivities within the 5-function approximation (Kamae et al.'s model A and Aharonian and Atoyan's 
formula for the inelastic cross section) with computations using differential cross section parameterizations. 



paper. It was found that a single expression is needed to repre- 
sent the total inclusive cross section, 

ln(£'„ - m„) 0.3 , 

crAEp) = (0.007 + 0.1-^,^ f + ^)-'mb, (9) 

(Ep - nip) (Ep - rripY 

where rest masses and energies must be given in units of GeV. 
These two differential cross section parameterizations pro- 
posed by Blattnig et al. are not deprived of problems if extrap- 
olated to high energy. The parameterization of the Stephen and 
Badhwar's model grossly underpredicts, whereas the newest 
Blattnig et al. (equation|S) overpredicts, the highest energy pion 
yield. The y-ray photon yield that is output of the use of these 
two differential cross sections in Equations ( 1415 > is also shown, 
for an Earth like spectrum, in FigureFTTI 

However, the inclusive total inelastic cross section (|9jl 
seems to work well at energies higher than 50 GeV, what we 
show in Figure ^] together with a compilation of experimen- 



tal data (Dermer 1986b). We also show in the same Figure the 
results for the inclusive total cross section from model A of 
Kamae et al. (2005), obtained from his figure 5. Indeed, Kamae 
et al.'s model produces a slightly higher cross section, although 
both correlate reasonably well with experimental data, at least 
up to 3 TeV^ It is the differential cross section parameterization 
(given by Equation|8|l the one that looks suspicious at such high 
energies. Figure^] right panel, shows the emissivities as com- 
puted in the different approaches multiplied by E^^^. As ex- 
pected, Stephen and Badhwar's parameterization falls quickly 
at high energy whereas both Kamae et al.'s and Aharonian and 
Atoyan's cross sections secure that the y-rays emitted maintain 



This figure enlarge the comparison of Blattnig et al. (2000) inclu- 
sive cross section with experimental data (see their figure 4), where 
only three low-energy data points from Whitmore (1974) were con- 
sidered. 
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an spectrum close to that of the proton primaries. For y-rays 
above few TeV, i.e., y-rays mostly generated by protons above 
few tens of TeV, Blattnig et al. differential cross section param- 
eterization makes the y-ray emitted spectrum much harder than 
the proton spectrum that produced them. This signals that a di- 
rect extrapolation of Eq. for computing photon emissivity 
above TeV with Eq. (|5} induces overpredictions of fluxes. 

Table 0] presents the results for the integrated emissivity, 
Q{E) dE, with Q(E) being the different curves of Figure 
lltl To obtain integrated fluxes from a source of volume V at 
a distance D one has to multiply by the constant V/ {4nD^), so 
that the difference in integrated emissivities indeed represent 
those among integrated fluxes. As Table0]shows -disregarding 
those coming from the Blattnig et al.'s parameterization of 
Stephen and Badwhar's results which are quoted here just for 
completeness- above 100 MeV, differences are less than a fac- 
tor of 1.5, which most likely is washed away by other uncer- 
tainties in any given model. But above 300 GeV, difference are 
larger and a conservative choice is in order. 

If interested in the GLAST-domain (say, £ > 100 MeV, 
£■ < 50 GeV) predictions, the most conservative choice seems 
to be the use of the Blattnig et al. new differential cross sec- 
tion parameterization (Eq. |8}, with no other approximation, in 
Equations and (|5}. This choice, while not taking into ac- 
count non-diffractive processes will possibly slightly under- 
predict the integrated flux (as shown in Table 0}. Up to this 
moment, there is no public parameterization of the differen- 
tial cross section including diffractive effects, but one is to be 
presented soon (T. Kamae, private communication). By using 
Blattnig et al. approach, there is no ^-function approximation 
involved nor an ad-hoc histogram of particle numbers as pro- 
posed in the treatment of Kamae et al. (2005). One has an an- 
alytical expression that can be directly used in the numerical 
estimates of Eq. (|5}. However, the price to pay is that this form 
of computation cannot be considered reliable at higher energies 
and should not be used. 

For the lACTs-domain (E > 100 GeV) the safest and also 
computationally-preferable choice appears to be to take either 
Kamae et al.'s model A, or even the simpler Aharonian and 
Atoyan's expression (Eq. and a 5-function approximation. 
This approach would probably be slightly underestimating the 
integrated flux at such high energies. All in all, assuming ei- 
ther Kamae et al.'s or Aharonian and Atoyan's expression for 
all energies does not import substantial differences, as Table 0] 
shows, and is computationally preferable. 

Finally, in Figure [O] we compare the inclusive cross sec- 
tions for charged pions with experimental data up to about 1 
TeV, which are again found in agreement with experimental 
data (except for a bunch of data points at low proton energies, 
in the case of positive charged pions). In any case, for situations 
where the density of cosmic rays or of target nuclei or both are 
high, neutral pion decay is expected to be the dominant process 
above 100 MeV, so that possible uncertainties in parameteriza- 
tions of charged pions cross sections are not expected to play a 
significant role in the prediction of fluxes. 
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Table 4. Integrated emissivities for an Earth-like spectrum. Values are in units of photons cm"^ s ' . 



Parameterization/ Approx. 


E > 100 Mev 


E > 100 GeV 


E > 315 GeV 


Blattnig et al. 


3.2E-25 


2.2E-28 


4.7E-29 


Kamae et al. 


4.8E-25 
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2.6E-30 


Aharonian 


4.0E-25 


2.2E-29 


3.1E-30 


Stephen and Badwhar (from Blattnig et al.) 


2.2E-25 


1.8E-31 


1.9E-33 




Fig. 13. Comparing inclusive cross sections for charged pions. Solid curves are obtained from Equations (30) and (3 1) of Blattnig 
et al.'s work (2000b), and experimental compilation is from Dermer (1986b). 
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